Preprint typeset using lAT^^X style emulateapj v. 04/03/99 



FORMING THE FIRST STARS IN THE UNIVERSE: THE FRAGMENTATION OF 

PRIMORDIAL GAS 

VoLKER Bromm, Paolo S. Coppi, and Richard B. Larson 

Department of Astronomy, Yale University, New Haven, CT 06520-8101; 
volker@astro.yale.edu, coppi@astro.yale.edu, Iarson@astro.yale.edu 






o 



> 

(N 

O 

OS- 

'^' 
Oh; 

6- 

Cd ■ 



% 



ABSTRACT 

In order to constrain the initial mass function (IMF) of the first generation of stars (Population III), 
we investigate the fragmentation properties of metal-free gas in the context of a hierarchical model of 
structure formation. We investigate the evolution of an isolated 3 cr peak of mass 2 x IO^Mq which 
collapses at Zcoii — 30 using Smoothed Particle Hydrodynamics. We find that the gas dissipatively 
settles into a rotationally supported disk which has a very filamentary morphology. The gas in these 
filaments is Jeans unstable with Mj ~ 10'^ Mq. Fragmentation leads to the formation of high density 
(n > 10^ cm^'^) clumps which subsequently grow in mass by accreting surrounding gas and by merging 
with other clumps up to masses of '^ W^Mq. This suggests that the very first stars were rather massive. 
We explore the complex dynamics of the merging and tidal disruption of these clumps by following their 
evolution over a few dynamical times. 

Subject headings: cosmology: theory — early universe — galaxies: formation — hydrodynamics 



1. INTRODUCTION 

Little is known about the history of the universe at red- 
shifts z ~ 1000 — 5, corresponding to 10^ — 10^ years after 
the Big Bang, and only recently have astrophysicists be- 
gun to seriously investigate this crucial post-recombination 
era. At z ~ 1000, the photons of the cosmic microwave 
background (CMB) shifted into the infrared, and the uni- 
verse entered what has been termed the "Dark Ages" (Rees 
1999). We know that the universe was reionized again be- 
fore a redshift of ^ 5 from the absence of Gunn-Peterson 
absorption in the spectra of high redshift quasars. One 
of the grand challenges in modern cosmology is to deter- 
mine when the universe was lit up again by the first lu- 
minous objects (the so-called Population III) which were 
responsible for reionizing and reheating the intergalactic 
medium (Haiman & Loeb 1997; Ferrara 1998). The very 
first generation of stars must have formed out of probably 
unmagnetized, pure H/He gas, since heavy elements can 
only be produced in the interior of stars. These character- 
istics render the primordial star formation problem very 
different from the present-day case, and lead to a signif- 
icant simplification of the relevant physics (Larson 1998; 
Loeb 1998). 

To determine the characteristic mass scale of the Popu- 
lation HI stars, and to constrain their initial mass function 
(IMF), one has to study the collapse and fragmentation of 
metal-free gas. This problem has been addressed by a 
number of authors (e.g., Yoneyama 1972; Hutchins 1976; 
Carlberg 1981; Kashlinsky & Rees 1983; Palla, Salpeter, 
& Stabler 1983; Silk 1983; Haiman, Thoul, & Loeb 1996; 
Uehara et al. 1996; Tegmark et al. 1997; Omukai & 
Nishi 1998; Nakamura & Umemura 1999). Recently, three- 
dimensional cosmological simulations have reached sufh- 
cient resolution on very small scales to approach the Pop- 
ulation HI star formation problem (Anninos & Norman 
1996; Ostriker & Gnedin 1996; Abel et al. 1998a). 

Complementary to these last studies, we explore the 
fragmentation of primordial gas under a variety of initial 



conditions. Our method allows us to follow the evolution 
of the fragments for a few dynamical times, and to include 
the physics of competitive accretion. These capabilities 
are crucial for the investigation of the resulting mass spec- 
trum. A more detailed exposition of our code, together 
with further results of our exploratory survey, will be pre- 
sented in a forthcoming publication (Bromm, Coppi, & 
Larson 1999). In this Letter, we focus on one represen- 
tative experiment which addresses the key question: How 
does the fragmentation of the first collapsing gas clouds 
proceed in a hierarchical (bottom up) scenario of structure 
formation ? 

2. NUMERICAL METHOD 

Our code is based on a version of TREESPH (Hernquist 
& Katz 1989) which combines the Smoothed Particle Hy- 
drodynamics (SPH) method (e.g., Monaghan 1992) with 
a hierarchical (tree) gravity solver. To study primordial 
gas, we have made a number of additions. Most impor- 
tantly, radiative cooling due to hydrogen molecules has 
been taken into account. In the absence of metals, H2 is 
the main coolant below ~ 10^ K, the typical temperature 
range in collapsing Population III objects. We implement 
the H2 cooling function given by Galli & Palla (1998). At 
temperatures approaching 10* K, cooling due to lines of 
atomic hydrogen dominates (Katz & Gunn 1991). The 
efficiency of H2 cooling is very sensitive to the H2 abun- 
dance. Therefore, it is necessary to compute the nonequi- 
librium evolution of the primordial chemistry. We take 
into account the six species H, H"*", H~, H2, HJ, and e~. 
Helium is neglected, since the temperatures in our appli- 
cation are low enough (< 10* K) to render the He species 
almost inert. Reaction rates are taken from Haiman et 
al. (1996). In addition, we have included 3-body reactions 
which become important at high density (n > 10^ cm~^), 
and can convert the gas almost completely into molecular 
form (Palla et al. 1983). The chemistry network is solved 
by the approximate backwards differencing formula (BDF) 
method (Anninos et al. 1997). 



We have devised an algorithm to merge SPH particles in 
high density regions to overcome the otherwise prohibitive 
timestep limitation, as enforced by the Courant stability 
criterion (Bate, Bonnell, & Price 1995). To follow the sim- 
ulation for a few dynamical times, we allow SPH particles 
to merge into more massive ones, provided they exceed 
a pre-determined density threshold, typically 10^ — 10^° 
cm^'^. More details of the merging algorithm will be given 
in Bromm et al. (1999). 

3. THE SIMULATION 
3.1. Initial Conditions 

Within a hierarchical cosmogony, the very first stars 
are expected to form out of 3 — 4cr peaks in the random 
field of primordial density fluctuations. The early (linear) 
evolution of such a peak, assumed to be an isolated and 
roughly spherical overdensity, can be described by the top- 
hat model (e.g., Padmanabhan 1993). We use the top-hat 
approximation to specify the initial dark matter (DM) con- 
figuration, where we choose the background universe to be 
described by the critical Einstein-de Sitter solution with 
density parameters: floM = 0.95, fls ~ 0.05, and Hub- 
ble constant h = iJo/(100 km s^^ Mpc"i)=0.5. In this 
paper, we investigate the fate of a 3a peak of total mass 
2 X 1O^M0, corresponding to 10^ Mq in baryons which is 
close to the cosmological Jeans mass (Couchman & Rees 
1986). On this mass scale, the standard Cold Dark Mat- 
ter (CDM) scenario predicts a present-day r.m.s. overden- 
sity of cro(Af) — 16, with a normalization erg = 1 on the 
8/i~^Mpc scale. Then, one can estimate the redshift of col- 
lapse (or virialization) from: 1+Zcoii = 3a-o{M)/l.Q9, lead- 
ing to Zcoii — 30. Our simulation is initialized at z^ = 100, 
by performing the following steps. 

The collisionless DM particles are placed on a cubi- 
cal Cartesian grid, and are then perturbed according to 
a given power spectrum P{k) = Ak", by applying the 
Zel'dovich approximation (Zel'dovich 1970) which also 
allows to self-consistently assign initial velocities. The 
power-law index is set to n = —3 which is the asymptotic 
small-scale behavior in the standard CDM model (Pee- 
bles 1993). The amplitude A is adjusted, so that the fun- 
damental mode (i.e., the smallest contributing wavenum- 
ber kmin), has an initial mean square overdensity which 
grows to P{kmm) — 1 at collapse. Next, particles within 
a (proper) radius of Ri — 150 pc are selected for the 
simulation. The resulting number of DM particles is 
Ndm = 14123. Finally, the particles are set into rigid 
rotation and are endowed with a uniform Hubble expan- 
sion (see also Katz 1991). Angular momentum is added by 
assuming a spin-parameter A — L\E\^^'^/{GM^/^) ~ 0.05, 
where L, E, and M are the total angular momentum, en- 
ergy, and mass, respectively. 

The colhsional SPH particles {Nsph = 16384 in this 
simulation) are randomly placed to approximate a uni- 
form initial density. The random sampling inevitably in- 
troduces shot-noise. The DM particles were set up on a 
regular grid specifically to avoid the unphysical shot-noise 
distribution which could mask the desired physical power 
spectrum. For the gas, however, the presence of the shot 
noise is not a big problem, since the gas mass is initially 
smaller than the Jeans mass. Therefore, sound waves will 
eflficiently wipe out all initial density disturbances. The 



SPH particles are endowed with the same Hubble expan- 
sion and rigid rotation as the DM ones. Since z < 100, 
heating, cooling, and photoreactions due to the CMB can 
be safely neglected here. For the chemical abundances 
and the gas temperature, we adopt the initial values (see 
Tegmark et al. 1997): x^ = 4.6 x 10"'', /^ = lO"'', and 
Tgas,i — 200 K, where Xi, fi are the initial free electron and 
H2 abundances, respectively. The use of a more realistic 
rate for the photodissociation of H2 predicts a lower ini- 
tial H2 abundance, typically fi ~ 10~®. The evolution of 
a simulation initialized with such a lower value, however, 
quickly converges to the case presented in this paper. 

3.2. Results 

3.2.1. Free-fall Phase 

From turnaround (z = 50) to the moment of collapse 
(z ~ 30), both the DM and the gas are freely falling. In re- 
sponse to the initially imprinted P{k) oc fc^'^ noise, the DM 
develops marked substructure. This happens in a hierar- 
chical fashion, where smaller clumps form first, and rapidly 
merge into larger aggregations. The gas initially falls to- 
gether with the DM, heating up adiabatically, T oc n?^^ 
for 7 = 5/3, until the temperature reaches the virial value 
of Tmr ^ 5000 K which is determined by the depth of 
the DM potential well. The virialized gas is able to cool 
efficiently, and the temperature decreases again with con- 
tinuing compression down to a few 100 K. At first, the 
gas is not able to condense onto the shallow DM potential 
wells, and its mass distribution remains smooth. In the 
course of the collapse, however, the baryonic Jeans mass 
decreases enough to eventually allow the gas to fall into the 
deepening DM troughs. This gravitational "head-start" is 
important because it determines where fragmentation oc- 
curs first, resulting in the most massive gas clumps. 

Briefly after collapse, all the substructure in the DM 
has been wiped out in the process of violent relaxation 
(Lyndcn-Bcll 1967) which is very effective in quickly es- 
tablishing virial equilibrium on the dynamical timescale 
of the DM halo (~ 10^ years). Before having been wiped 
out, however, the DM substructure has left its imprint on 
the gas, determining the morphology of the incipient disk. 
The DM reaches an equilibrium configuration with a core 
radius of ^ 10 pc and an approximately isothermal density 
profile [p oc r^^) further out. 

The central gas disk is horizontally supported by rota- 
tion. In the presence of a DM halo, one expects a contrac- 
tion by a factor of 1/A ~ 20 for rotational support. The di- 
mension of the disk, as shown in Figure 2, is in good accor- 
dance with this prediction. Vertically, the disk is roughly 
pressure supported with a scale height of iJ = c^/(GE) ~ 
2 pc, where c^ ~ 2 km/s is the sound speed, and S the 
gas surface density. The resulting disk has a very filamen- 
tary and knotty structure with typical densities of n '-~^ 10 
cni^^ and temperatures of T ^^ 500 K, corresponding to 
a Jeans mass of Mj ~ IO^Mq. We next discuss the frag- 
mentation of these filamentary features. 

3.2.2. Disk-like Phase 

The gas which makes up the filaments is Jeans unsta- 
ble since tff < t sound, where the free-fall time is tff ~ 
{Gmiin)~^'^ , and the sound-crossing time t sound — R/cs- 
i? ~ 15 pc is the radius of the disk. In Figure 1, the free- 



fall time is compared to the sound-crossing time and to the 
cooling time, tcooi — nkBT/{fn^L). Here, / is the frac- 
tional H2 abundance, and L the cooling rate (in erg s~^ 
cm^). For densities n > 10"^ cm~^, the gas becomes Jeans 
unstable. The onset of instability roughly coincides with 
the condition that tcooi — iff. The value n ~ 10"^ cm~'^ is 
also the critical density for the rotational- vibrational levels 
of H2 being populated according to local thermodynamic 
equilibrium (LTE). 

The first regions to undergo runaway collapse lead to 
the formation of dense (n > 10^ cm~^) gas clumps of 
mass ~ 10^ Mq. The dynamical state of a typical clump 
can be described by the ratios Eth/\Egrav\ — 0.5, and 
Erot/\Egrav\ - 0.1, whcrc Efh, Erot, Egrav are the ther- 
mal, rotational, and gravitational energies of the clump. 
The further evolution of the disk consists of a complex 
history of new clumps forming, and old ones continually 
accreting gas and occasionally merging with other clumps. 
The outcome of this evolution is shown in Figure 2, corre- 
sponding to the situation at z = 28. The disk- like struc- 
ture has fragmented into 13 clumps, ranging in mass from 
220 to ~ IO^'Mq. At this stage, --50% of the total gas mass 
has been incorporated into high density clumps. This frac- 
tion is unrealistically high due to the homologous nature 
of the top-hat collapse where all particles reach the cen- 
ter at the same time. Abel et al. (1998a), on the other 
hand, estimate this fraction to be '^8%. This value is 
more realistic, since these authors consider an overdensity 
which is initially already centrally concentrated, leading 
to a collapse which is spread out in time. Since the larger 
clumps tend to retain their individuality, the final evolu- 
tionary state of the disk is expected to be a small cluster 
of isolated clumps. 

Figure 3 summarizes the thermodynamic state of the gas 
together with the chemical abundances. When a region 
undergoes runaway collapse, particles are drawn from the 
quasi-stationary "reservoir" at n '^ 10* cm^'^ which corre- 
sponds to the filamentary gas, and swiftly move to higher 
densities at roughly constant temperature (Fig. 3c). At 
a density of n = lO^cm^'^, SPH particles are merged into 
more massive ones, as described in §2. At higher densities, 
the Jeans mass would decrease below the resolution limit 
{Mres ^ 2OOM0; Bate & Burkert 1997) of our simulation 
(Fig. 3d). In order to reliably follow the evolution of the 
gas to even higher densities, a larger number of SPH par- 
ticles would be needed. In this case, one would expect the 
gas density to increase up to the point where the now fully 
molecular gas becomes opaque, and subsequently behaves 
adiabatically again. The present study, however, has the 
purpose of following the evolution of many gas clumps for 
a few dynamical times. Therefore, merging at some point 



is inevitable. 

4. SUMMARY AND CONCLUSIONS 

We have discussed the evolution of a primordial star 
forming region in the context of the standard CDM model 
for structure formation. In the early stages of the collapse, 
the morphology of the DM substructure imprints a pattern 
onto the initially completely smooth gas which determines 
where gas fragmentation is to occur first. Therefore, the 
initial power spectrum of the DM fluctuations might be 
an important ingredient for primordial star formation. Si- 
multaneously with the virialization of the DM component, 
the gas dissipatively settles down into a compact central 
disk. This disk has a very filamentary structure, and is 
Jeans unstable which leads to fragmentation into clumps 
of initially ~ 10^ — IO^Mq. These clumps grow in mass 
by both accreting surrounding gas, and by merging with 
other high-density clumps. At the end of this process, a 
small number ('^ 10) of these clumps have formed. 

Population HI stars, then, might have been quite mas- 
sive, perhaps even 'very massive' (Mchar > IOOMq), if one 
assumes that a Jeans mass clump does not undergo signif- 
icant further fragmentation. Although subsequent frag- 
mentation cannot be excluded, this is expected to lead 
to the formation of a binary or a small multiple system, 
if experience from present-day star formation offers any 
guidance here (Larson 1998). Our results agree well with 
those of Abel, Bryan, & Norman (1998b) who have ap- 
plied the Adaptive Mesh Refinement (AMR) method to 
the problem. These authors find a dense core of order a 
few times IO^Mq, which continues to accrete mass and 
shows no sign of further subfragmentation. Clearly, high 
resolution studies of the further clump evolution are neces- 
sary, including a treatment of opacity effects (see Omukai 
& Nishi 1998). We are presently engaged in doing so and 
will present our findings in a future publication. 

The effects of cooling due to HD molecules (Galli & Palla 
1998), and of the presence of the CMB on the H2 level pop- 
ulations (Capuzzo-Dolcetta, Di Fazio, & Palla 1991), were 
not included in our work, and both could affect the esti- 
mate of the Jeans mass. Test calculations indicate that the 
effect of HD cooling may not significantly change our re- 
sults, but a more complete treatment will be incorporated 
in subsequent work. 
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Palla for detailed comments. Support from the NASA 
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Fig. 1. — Important timescales at z = 28. (a) Ratio of cooling and free-fall tiineseales vs. gas density. At low density, cooling is not 
efficient and contraction proceeds adiabatically. For higher densities, characteristic of the gas after virialization, cooling is very efficient up 
to densities n ~ 10^ cm~^. The evolution from n ~ 10^ cm~^ to n ~ 10* cni""^ proceeds such that t^ool — iff- (b) Ratio of sound-crossing 
and free-fall timescales vs. gas density. For tff < tgound, the gas is Jeans unstable. Gravitational instabilty sets in at n ^ 10 cm^ which 
corresponds to the compressed gas in the disk-like configuration. The onset of instability roughly coincides with the condition that t^^^i ~ tff 
in panel (a). 
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Fig. 2. — Situation at z = 28. Shown is the ongoing fragmentation of the gas disk in the center of the DM halo. The (proper) length of 
the box is 30 pc. Left Panels: Face-on view. Right Panels: Edge-on view, (a)-(b) Gas component. The gas which had been assembled in 
a disk-like configuration, has undergone further fragmentation into roughly spherical high-density (n > 10* cm~^) clumps, (c)-(d) Clump 
distribution. The four increasing dot sizes denote increasing mass scale: 10^ — 10^ Mq, 10^ — 5 X 10^ Mq, 5 X 10"^ — 10* M©, and M > IO^Mq. 
By now, '^ 50 % of the total gas mass is incorporated into high-density clumps. The most massive clumps have already accreted all of the 
surrounding gas. This explains the paucity of unmerged gas particles close to the center of the disk in panel (a). 
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Fig. 3. — Gas properties at z = 30.5. (a) Free electron abundance vs. gas number density. At high density, the residual electrons recombine 



until the gas is almost completely neutral (at n -^ 10'' 
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abundance vs. gas number density. Due to the H channel, the 
(c) Gas temperature vs. number density. At low gas densities, the 



temperature rises due to adiabatic compression until it reaches the virial value. At higher densities, H2 line cooling drives the temperature 
down again until the gas settles into a quasi-equilibrium state at T ~ 500 K and n ~ 10 cm~ . (d) Jeans mass vs. number density. The Jeans 
mass reaches a value of Mj ~ 10^ Mq in the disk-like central feature, and approaches the resolution limit of the simulation, Mres ~ 200Mq, 
for densities close to the merging threshold of n = 10* cm~^. 



